% ***************************************************************************************
% m_calibration_1930s.m 
% ***************************************************************************************
clear all; 
close all;
format shortG;

% ==================== SET FUNDAMENTAL PARAMETERS ========================
load parameters.mat;                     % basic parameters 
load calib_productivity_amenity.mat;     % inverted fundamentals every 5 years (55,60,65,70,75)
load commuting_cost.mat;  

% =================== PARAMETERS ==========================================
alpha=alphaest;   
beta=betaest; 
kappa30s=Kappa_mat(:,:,1); 
K30s=kappa30s.^(-rho50/sigma50); 

% ========================= LOAD DATA for 1930 =========================================  
datafile = '../../data/processed data/quant/Hiroshima_DATA.csv';
CITYDATA = readtable(datafile);  
ID = CITYDATA.geocode;  
N=size(ID,1); 
Area = CITYDATA.area; 
GEOGRAPHY=[ID Area]; 
cbddist=CITYDATA.cbddist;

Rdata36=CITYDATA.pop36; Pop36=sum(Rdata36); 
Ldata38=CITYDATA.emp38; Emp38=sum(Ldata38); 

% ************** ADJUSTMENT FOR 1930S BY SCALE UP FOR POPULATION **********
R30s=Rdata36; 
L30s=Ldata38.*(Pop36/Emp38); 
if (sum(R30s) ~= sum(L30s)) == 1   % check total pop ; 
    disp('>>> CHECK TOTAL POPULATION AND EMPLOYMENT ! <<< ');   pause
end 

% ******************* GRAVITY FOR COMMUTING *******************************
X0=ones(N,1);
options0 = optimset('Display','iter','Algorithm','trust-region-dogleg','MaxFunEvals',50000000,...
          'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
profile on 
syst0 = @(X)fcomeval(X,kappa30s,R30s,L30s,gravity);
[x_fsolve, fval_fsolve]=fsolve(syst0, X0, options0);
w30s=abs(x_fsolve);
w30s=w30s./geomean(w30s);

% ****************** AVERAGE INCOME IN RESIDENCE **************************
lambdatemp=(repmat(w30s',N,1)./kappa30s).^gravity;
lambdacon=lambdatemp./repmat(sum(lambdatemp,2),1,N);  
lambda=lambdacon.*repmat(R30s./sum(R30s),1,N);        
v30s=lambdacon*w30s;                                  

% ****************** FLOOR SPACE MARKET CLEARING CONDITION ****************
floordemand30s = mu.*v30s.*R30s + (gammaH/gammaL).*w30s.*L30s; 
A30s = (w30s.^gammaL).*(psi.*floordemand30s./Area).^(1/(1+eta/gammaH));
A30s = A30s./geomean(A30s); 
Q30s = (A30s./(w30s.^gammaL)).^(1/gammaH); 

% ****************** RESIDENTIAL AMENITIES ********************************
options0=optimset('Display','iter','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
B0=ones(N,1); 
Rsyst = @(x)fReval(x,K30s,R30s,L30s,rho50,sigma50);
[x_fsolve, xfval_fsolve]=fsolve(Rsyst, B0, options0); 
B30s=abs(x_fsolve); 
B30s=B30s.*(Q30s.^mu);
B30s=B30s./geomean(B30s);

% ****************** DECOMPOSE FUNDAMENTALS *******************************
a30s=A30s.*((L30s./Area).^(-alpha)); 
b30s=B30s.*((R30s./Area).^(-beta)); 

% ****************** SAVE FUNDAMENTALS FOR 1930S **************************
fundamentalresult=array2table([w30s,Q30s,a30s,b30s,aa,bb,Rdata36, Ldata38, cbddist, Area, ID],'VariableNames', {'w30','q30','aa30',...
    'bb30', 'aa55','aa60','aa65','aa70','aa75','bb55','bb60','bb65','bb70','bb75', ...
    'pop36','emp38','cbddist', 'area', 'geocode'}); 

writetable(fundamentalresult,'../../output/quant/output_calib_fundamentals_1930s.csv');

% ************************************************************************
% ************************************************************************
% **************************  FUNCTIONS **********************************
% ************************************************************************
% ************************************************************************

% --------  FUNCTION FOR COMPUTING WAGE BASED ON GRAVITY STURCTURE --------
function[commuteeval] = fcomeval(X,CD,DR,DL,epsilon)
N = size(DR,1); 
X=abs(X);  
X=X./geomean(X);
lambdatemp=(repmat(X',N,1)./CD).^epsilon;
lambda=lambdatemp./repmat(sum(lambdatemp,2),1,N);
commutedif = DL-lambda'*DR; 
commuteeval = abs(commutedif);
end

% --------  FUNCTION FOR COMPUTING CHARACTERISTICS BASED ON GRAVITY STURCTURE --------
function[Reval] = fReval(B,K,R,L,rho,sigma)  
N = size(R,1); 
B=abs(B); 
B=B./geomean(B);
xtemp = K.*(repmat(B',N,1).^(rho/sigma)); 
xtemp = xtemp./repmat(sum(xtemp,2),1,N); 
Reval = R-sum(xtemp.*repmat(L,1,N),1)'; 
Reval = abs(Reval); 
end

